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ABSTRACT 

Significant  contributions  in  algorithms  and  codes  for 
automatic  grid  generation  were  made  in  the  following 
specific  areas:  automatic  generation  of  multidimensional 
Fortran  code  via  Symbolic  Manipulation  (Artificial 
Intelligence) ,  requiring  user  input  of  only  the  governing 
partial  differential  equations  or  the  governing  variational 
principles;  rigorous  techniques  for  verification  of  codes, 
algorithms,  and  discretization  methods;  a  new  and  effective 
formulation  for  variational  methods;  analysis  of 
mathematical  properties  (ellipticity,  convexity)  of  PDE  and 
variational  methods  including  volume,  orthogonality, 
smoothness,  and  other  controls;  identification  and 
explanation  of  previously  unreported  anamolies  in  common 
grid  generation  algorithms,  including  non-uniqueness  due  to 
solution  bifurcation  and  grid  folding;  identification  and 
explanation  of  previously  unreported  ambiguities  in  common 
discrete  evaluations  of  transformation  Jacobians;  use  of  a 
reference  grid  to  control  target  grid  properties; 
development  of  a  solution  adaptive  method  for  surface  grid 
generation;  development  of  a  hybrid  method  of  solution 
adaptive  grid  generation  which  properly  isolates  the 
contribution  of  the  adaptivity  functional;  application  of 
multigrid  methods  to  elliptic  grid  generation;  development 
of  an  algorithm  for  economically  generating  sub-grid  grid 
generation  discretizations  required  for  multigrid  solution 
techniques;  correct  interpretation  of  the  "smoothness" 
property  of  elliptic  grid  generation. 
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B.  STATEMENT  OF  WORK 

"The  contractor  will  develop,  verify  and  exercise 
Symbolic  Manipulation  codes  and  Fortran  codes  for 
solution-adaptive  grid  generation  for  two  dimensional  and 
three  dimensinoal  internal  and  external  flows.  The  methods 
and  results  produced  under  the  proposed  contract  will  be 
documented  in  open  literature  and  engineering  meetings,  and 
sample  codes  will  be  made  available  to  interested  users 
within  Government  agencies  and  their  contractors." 

C.  STATUS  OF  THE  RESEARCH,  and  D.  TECHNICAL  JOURNAL 
PUBLICATIONS . 


The  algorithms,  mathematical  analyses,  and  engineering 
applications  outlined  herein  were  developed  under  the 
subject  contract  and  its  predecessor  (F49620-82-C-0064 ,  same 
title) .  The  original  intention  of  AFOSR  was  to  simply 
extend  the  earlier  contract;  this  failed  only  because  of  a 
backlog  at  the  AFOSR  contracts  office  at  the  time. 
Consequently,  the  division  of  results  between  the  two 
contracts  is  somewhat  vague  scientifically.  Reference 
should  also  be  made  to  the  earlier  work  reported  in  the 
Ecodynamics  Final  Scientific  Report  on  the  previous  contract 
dated  10  May  1986, 


In  addition  to  journal  publications,  two  books 
eventually  resulted  from  the  work  begun  under  sponsorship  of 
the  subject  contract  (as  well  as  collateral  funding) . 

The  first  book  is  "Mathematical  Aspects  of  Numerical 
Grid  Generation",  No.  8  in  the  SIAM  Series  "Frontiers  in 
Applied  Mathematics",  edited  by  Dr.  Jose'  E.  Castillo. 

Dr.  Castillo  was  partially  supported  by  the  subject 
contract  as  an  Ecodynamics  employee  while  he  was 
successfully  pursuing  his  Ph.D.  in  Mathematics.  The 
Principal  Investigator,  Dr.  P.  J.  Roache,  and  Dr.  S. 
Steinberg  were  his  co-advisors.  Five  of  the  ten  chapters  of 
this  book  were  written  by  Ecodynamics  personnel . 


The  second  book  is  "Fundamentals  of  Grid  Generation"  by 
P.  M.  Knupp  and  S.  Steinberg.  Still  in  manuscript  form, 
it  is  expected  to  be  published  in  1993,  probably  by  SIAM. 


The  most  economical  and  meaningful  method  of 
summarizing  the  results  of  the  contract  work  is  to  simply 
duplicate  herein  the  title  page,  including  the  technical 
abstract,  for  the  publications  supported  by  the  subject 
contract.  These  follow. 
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Abstract 

This  paper  reports  on  three  aspects  of  grid  genera¬ 
tion.  First,  two  behavioral  errors  in  elliptic  gird  gen¬ 
eration  that  are  of  critical  importance  in  choosing  a 
grid  generation  method  are  discussed.  Next,  a 
geometric  continuation  method  for  avoiding  difficulties 
sometimes  encountered  in  difficult  geometries  is 
described.  Finally,  a  significant  extension  of  the  exist¬ 
ing  variational  grid  generation  methods  is  described. 

Background 

We  present  here  a  brief  background  of  elliptic  grid 
generation  relevant  to  the  present  work.  Winslow 
(Ref.  1)  used  homogeneous  Laplace  equations  written 
in  the  physical  plane  for  grid  generation.  These  equa¬ 
tions  are  transformed  to  the  logical  plane  where  they 
are  solved  numerically.  The  concept  behind  this 
approach  was  based  on  the  utilization  of  the  well- 
known  maximum  principle  for  linear  elliptic  equations. 
The  maximum  principle  for  the  Laplace  equation 
guaranteed  that  the  method  would  always  produce  a 
minimally  acceptable  grid,  i.e.  one  for  which  the  grid 
lines  would  not  leave  the  region.  Solving  the  equations 
in  the  transformed,  or  logical  plane,  is  advantageous 
due  to  a  simplification  of  the  boundary  conditions. 
However,  this  occurs  at  the  expense  of  a  nonlinear  cou¬ 
pling  of  the  governing  equation.  Winslow  used  a 
discretization  on  a  triangular  grid. 

Hirt,  .\msden.  and  Cook  (Ref.  2)  used  Laplace 
equations  written  in  the  logical  plane  itself.  Here,  the 
maximum  principle  does  not  help,  so  the  continuum 
formulation  has  no  guarantee  of  a  minimally  accept¬ 
able  grid.  However,  the  method  is  linear  and  the 
Laplace  equations  are  uncoupled.  The  method  was  able 
to  produce  usable  solutions  in  some  cases. 

Thompson,  Thames,  and  .Mastin  (TTM  method. 
Ref.  3)  extended  the  Winslow  method  by  introduction 
of  nonhomogeneous  terms  in  the  physical  plane,  fol¬ 
lowed  by  transformation  to  the  logical  plane.  This 
resulted  in  a  much  more  complicated  nonlinear  equa¬ 
tions,  but  the  right-hand  side  or  nonhomogeneous 
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terms  allows  control  of  the  grid  at  interior  points.  The 
control  achieved  by  tbe  nonhomogeneous  terms  is  very 
indirect;  several  papers  have  been  published  on  the 
forms  of  the  nonhomogeneous  terms,  which  are  strong 
and  typically  greatly  slow  the  iterative  convergence 
rate.  See,  for  example,  Sorenson  and  Steger  (Ref.  4), 
Middlecoff  and  Thomas  (Ref.  5),  Ghia  et  al  (Ref.  8) 
and  Brackbill  and  Saltzman  (Ref.  7).  We  note  that  the 
introduction  of  the  nonhomogeneous  terms  means  that 
the  maximum  principle  no  longer  applies,  especially  if 
the  sign  of  the  nonhomogeneous  term  changes  within 
the  domain. 

Brackbill  and  Saltzman  (Ref.  7)  approached  the 
grid  generation  problem  using  a  variational  formula¬ 
tion  of  the  three  properties  of  smoothness,  cell  volume, 
and  orthogonality,  written  in  the  physical  plane  and 
then  transformed  to  the  logical  plane.  The  variational 
principal  for  smoothness  alone,  properly  defined,  leads 
to  the  Winslow  equations  (or  the  homogeneous 
Thompson-Thames-Mastin  equations).  Weighting  fac¬ 
tors  are  introduced  to  allow  the  user  to  balance  the 
importance  of  these  three  properties.  The  particular 
weighting  of  orthogonality  used  in  Ref.  7  b  somewhat 
arbitrary,  but  simplifies  the  mathematics  of  the 
transformation. 

Behavioral  Errors  in  Grid  Generation 

We  have  encountered  some  -.lementary  aspects  of 
elliptic  grid  generation  which  bear  on  the  selection  of 
an  existing  method  and  on  the  direction  of  the 
research  for  our  new  meth  ds. 

These  errors  fall  urder  the  category  of  “behavioral 
errors"  (Ref.  8),  i.e.  qualitative  differences  between  the 
dis'-retized  equations  and  the  continuum  equations. 
Thb  genre  of  classification  is  well-known  in  CFD  and 
includes  (ma_.)  conservation  errors,  vorticity  and 
enstrophy  enservation  errors,  phase  errors,  artificial 
vbcosiiy  errors,  positivity  errors  (e.g.,  local  appearance 
of  negative  values  for  positive  definite  qualities  such  as 
temperature  and  density).  Other  branches  of  computa¬ 
tional  physics,  e.g.,  electromagnetics,  experience  analo¬ 
gous  errors. 

The  first  behavioral  error  b  in  tbe  evaluation  of  cell 
volume  and  grid  folding  using  a  dberetized  Jacobian. 
Let  the  mapping  from  the  logical  plane  to  the  physical 
plane  be  given  by  *  =  .?((,7),  9  =  ir({,q).  In  2D,  the 
sign  of  the  Jacobian, 
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Abstract 

The  goal  of  the  computational  effort  described 
herein  is  to  develop  computer  codes  for  rapidly  and 
accurately  modeling  the  electric  fields  in  the 
cavity  of  lasers  and  switches.  The  designer  is  able 
to  interactively  perturb  the  laser  operating  para¬ 
meters  and/or  electrode  geometry,  and  ouickly  od- 
tain  new  solutions.  The  cooes  use  automatic 
generation  of  solution-adaptive  bounoary-fitted 
coordinate  systems,  and  solve  tv/o-  anc  tnree- 
dimensional  proolems  in  Doth  steaoy-state  and  time- 
depenoent  modes. 

Introduction  and  Overview 

The  design  of  electrodes  for  lasers  and  switches 
is  well-defined  only  for  unrealistically  idealized 
conditions.  The  frequently  used  Rogowski  electrode 
shapes  are  "optimal"  only  in  the  sense  of  producing 
an  enhancement  factor  of  unity,  i.e.,  the  electric 
field  strength  is  no  where  greater  than  the  nominal 
value.  More  importantly,  the  solution  is  based  on 
vacuum  conditions  and  is  not  a  complete  specifica¬ 
tion,  i.e.,  the  Rogowski  shape  is  not  closed,  and 
must  be  completed  by  some  (usually  arbitrary) 
closure  such  as  blending  with  a  radius,  etc.  The 
same  is  true  of  the  Chang  electrodes. 

The  computer  codes  described  herein  address  the 
realistic  electrode  design  problem,  including  non¬ 
vacuum  operation  and  complete  electrode  specifica¬ 
tion  with  "packaging"  constraints  of  overall  size. 
Using  efficient  finite  difference  methods  in 
boundary-fitted  coordinates,  the  ELF  (ELectric 
Field)  code  makes  it  practical  to  design  the  elec¬ 
trode  geometry  and  laser  operating  parameters 
during  user-interactive  sessions  on  a  VAX  computer. 

A  single  code  is  used  for  all  2-0  calculations, 
both  steady  state  and  time  dependent.  Options  are 
available  for  either  the  planar,  axisymmetric,  or 
radial  electrode  geometries.  Boundary  conditions 
and  boundary  snapes  may  be  time  dependent;  in  par¬ 
ticular,  an  external  circuit  equation  is  provided 
so  that  electrode  potential  may  be  calculated  as 
part  of  the  solution,  dependent  on  tne  integrated 
current  throuqn  the  cavity,  rather  than  being 
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Weapons  Laooratory,  the  U.S.  Air  Force  Office  of 
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specified  a  priori  .  The  geometry  and  conductivity 
calculations  are  modularized  so  that  they  may  be 
readily  modified  by  the  user. 

Automatic  grid  generation  is  performed  interac¬ 
tively  using  elliptic  generating  equation  techniques. 
As  an  option,  a  solution  adaptive  grid  generation 
technioue  is  used  to  adapt  tne  grid  to  tne  solution 
(either  in  tne  steady  state  solution,  or  within  an 
intra-time-step  iteration  for  a  time-dependent 
proDlem)  so  as  to  increase  the  resolution  of  tne 
maximum  electric  field  strength  (always  an  impor¬ 
tant  design  parameter)  and  the  accuracy. 

The  code  accounts  for  externally  controlled  or 
self-sustained  glow  discharges  or  other  plasmas, 
such  as  arcs,  by  modifying  the  nonlinear  conducti¬ 
vity.  Various  empirical  formulations  for  (steady) 
electrical  conductivity  have  been  evaluated  using 
the  code.  While  somewhat  useful,  this  approach  has 
now  been  abandoned.  A  more  fundamental  approach 
has  been  adopted,  that  of  solving  for  the  conducti¬ 
vity  by  time  integration  of  the  ordinary  differen¬ 
tial  equations  for  electron  numoer  density  equation 
at  each  mesh  point  in  the  2D  or  30  grid,  coupled 
nonlinearly  to  the  local  E-field,  The  electron 
drift  velocities,  etc.,  are  obtained  by  table  look¬ 
up  of  Boltzmann  code  solutions  performed  beforehand 
(i.e.,  noninteractively)  for  the  particular  gas 
mixture  used.  Code  applications  have  included 
pulsed  electric  CO2  lasers  and  a  xenon  flashlamp. 

For  xenon  flashlamp  calculations  and  for  streamer 
calculations,  the  temperature  is  also  solved  at 
each  node  point  by  time  integration  of  an  energy 
equation.  Calculations  have  shown  some  insight  into 
streamer  formation  in  plasma  discharges,  the  un¬ 
steady  development  of  a  self-sustained  glow  dis¬ 
charge,  and  lensing  effects  due  to  nonuniformities 
in  external  ionization  sources. 

For  many  cases  studied,  we  find  the  electric 
field  solutions  to  differ  significantly  from  vacuum 
calculations,  indicating  that  the  commonly  used 
Rogowski  solutions  and  Chang  solutions  for  the 
electrode  shapes  are  far  from  optimal  for  important 
classes  of  problems.  The  true  optimal  geometry  is, 
in  fact,  a  strong  function  of  the  lase*-  physics  and 
of  the  operating  conditions  whenever  significant 
physics  are  involved  in  the  conductivity.  Also, 
different  lasers  may  have  different  optimality  cri¬ 
teria;  e.g.,  an  el ectron-oeam  laser  may  be  designed 
to  give  nearly  uniform  energy  deposition  in  the 
cavity,  whereas  a  self-sustaineo  CO2  laser  may  be 
designed  to  minimize  the  local  extrema  of  the  elec¬ 
tric  field  strength,  subject  to  external  packaging 
geometry  constraints,  so  as  to  minimize  arcing  and 
maximize  the  operating  voltage. 

Three-dimensional  calculations  are  done  in  a 
separate  code  and  are  used  only  to  design  the  roll¬ 
off  of  the  electrooes  in  the  third  dimension  so  as 
not  to  produce  a  locally  high  electric  field  due  to 
3-0  effects. 
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capabilities  of  the  sstnholic  mantpulation  code  MACs^MA  deseloped  at  MIT.  The  second  section  presents 
results  recentls  achiesed  h>  the  authors  iisini;  symbolic  manipulation  on  the  problem  of  three-dimensional  hody- 
fittcd  coordinates.  The  third  part  of  the  paper  surveys  previous  uses  of  symbolic  manipulation  in  computational 
lluid  d>  namics  and  related  areas,  and  potenttai  future  uses  of  the  code. 


Introduction 

HE  purpose  ot'  thi.s  paper  is  to  provide  a  brief  in¬ 
troduction  to  the  uses  of  symbolic  manipulation  syithin 
the  practice  of  disciplines  such  as  computational  lluid 
dynamics  (CFD).  The  first  section  is  a  brief  presentation  of 
the  powerful  capabilities  of  the  symbolic  manipulation  code 
M.ACSYMA.  The  next  section  presents  results  recently 
achieved  by  the  present  authors  using  symbolic  manipulation 
on  the  problem  of  three-dimensional  boundary-filled 
coordinates.  The  last  section  consists  of  a  short  survey  of 
previous  uses  of  symbolic  manipulation  in  CFD  and  related 
areas,  and  of  what  we  consider  to  be  promising  areas  for 
future  use. 

Symbolic  manipulation  performed  by  computer  is  an  area 
that  has  just  recently  begun  to  be  more  widely  appreciated  and 
has  tremendous  potential.  Symbolic  manipulations  are  not 
floating  point  calculattons.  but  symbolic  operattons  (e.g.,  the 
chain  rule  differenttation)  performed  by  computer  logic.  One 
of  the  earliest  symbolic  manipulation  languages  was  FOR- 
M.AC.  which  had  the  capability  of  doing  simple  dif¬ 
ferentiation  symbolically.  .ALTR.AN  is  a  more  powerful  and 
later  version  of  the  symbolic  manipulation  language  another 
version.  "Micromath."  is  available  on  the  better 
microcomputers.  Perhaps  the  most  powerful  of  the  symbolic 
manipulation  codes  is  .M.ACSYMA.  which  has  been  developed 
over  many  years  at  the  .Massachusetts  Institute  of  Tech¬ 
nology. 

M.ACSYMA  Capubilities 

The  symbolic  manipulation  code  .MACSI'M.A  has  been 
developed  through  an  extensive  team  effort.  Many  prac¬ 
titioners  of  computational  lluid  dynamics  have  some 
acuiiainiance.  at  least  second-hand,  with  this  code,  but  our 
impression  is  that  few  have  an  appreciation  for  its  powerful 
capabilities  or  those  of  the  general  field  of  symbolic 
manipulation.  These  capabilities  include  variable-precision 
arithmetic  and  algebraic  substitutions  and  simplifications, 
svmbolic  (rather  than  numeric!  equation  solving  for  scalars 
and  matrices,  truncated  Taylor  series  expansions,  power 
series,  differentiation  including  the  chain  rule  for  single-  and 
multivariate  functions,  very  powerful  definite  and  indefinite 
integration,  taking  limits  of  expressions,  ordinary  differential 
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equation  solving  iclosed-form  solutions),  and  many  other 
mathematical  operations  previously  thought  to  be  in  the  realm 
of  human  beings  rather  than  computers.  For  an  introduction 
to  the  subtect  of  symbolic  manipulation,  see  Refs.  1-3. 
Recently,  a  version  of  MACSYM.A.  called  VA.XIMA,  has 
become  available  for  the  V.A.X  family  of  minicomputers. 
Significantly,  MACSY.MA  also  has  the  capability  of  actually 
generating  a  Fortran  code  from  symbolic  mathematics  in  the 
usual  scientific  notation.  MACSYM.A  is  also  a  programming 
language  and  can  concatenate  alphanumeric  variables. 

The  second  author  has  developed  a  code  called  GEEVVHIZ 
that  demonstrates  many  of  the  impressive  capabilities  of 
M.ACSYMA.  .A  sampling  is  included  as  an  appendix  in  Ref. 
21. 

Three-Dimensional  Boundarv-Fiffed  Coordinates 
via  Symbolic  .Manipulation 

Background 

Coding  errors  plague  all  computational  work.  The  chance 
for  coding  error  increases  as  the  complexities  of  the  problems 
increase.  Computational  lluid  dynamics  codes  used  in 
aerodynamics  and  ballistics  calculations  are  commonly  5,000 
or  more  lines  of  Fortran  code  in  length.  Codes  used  in  nuclear 
reactor  .safety  work  have  as  many  as  65,000  lines  of  Fortran. 
.Although  containing  much  complex  physics,  these  codes  are 
actually  elementary  in  that  they  are  written  in  very  simple 
coordinate  systems  such  as  Cartesian,  cylindrical,  etc. 

Recently,  there  has  been  a  surge  of  activity  in  the  ap¬ 
plication  of  coordinate  transformation  techniques.  The 
proceedings  of  three  recent  conferences,  two  dedicated  en¬ 
tirely  to  this  problem*-  and  one  emphasizing  it.'’  show  a 
tremendous  amount  of  work  being  done  in  this  area.  We  have 
been  impressed  with  the  difficulty  of  code  verification  for 
transformed  grid  problems  and  with  the  complexity  of  ihree- 
dmensional  equations  written  in  general  nonorthogonal  grids. 
The  most  popular  and  most  generally  applicable  coordinate 
transformation  methods  are  boundary-fitted  coordinate 
systems,  which  are  generally  nonconformal  and  therefore  do 
not  maintain  the  form  of  the  original  equations  written  in 
Cartesian  or  other  elementary  coordinates. 

Overview 

We  have  recently  developed  and  validated  symbolic 
manipulation  codes  for  working  a  problem  of  considerable 
comple.xity — that  of  solving  general  second-order  elliptic 
partial  differential  equations  in  a  general  (nonorthogonal) 
.three-dimensional  boundary-fitted  coordinate  system.  We 
follow  the  approach  of  Thompson  and  his  colleagues’  *  and 
generate  the  coordinate  system  itself  by  solving  elliptic  partial 
differential  equations  in  the  transformed  plane.  We  address 
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The  problem  of  numerically  integrating  general  elliptic  differential  equations  in  irregular 
two  and  three  dimensional  regions  is  discussed.  The  method  used  numerically  computes  a 
transformation  of  the  given  region  into  a  rectangular  region.  The  numerical  coordinate 
transformation  is  determined  by  requiring  that  the  components  of  the  transformation  satisfy 
inhomogeneous  Laplace  or  more  general  equations.  The  transformation  is  then  used  to 
transform  the  differential  equation  and  the  boundary  conditions  to  the  rectangular  region.  The 
boundary  value  problem  in  the  rectangular  region  is  integrated  using  one  of  the  standard 
methods  for  general  elliptic  equations.  The  use  of  the  existing  software  reduces  the  problem  to 
analytically  transforming  the  given  differential  equation  and  the  Laplacian  to  general  coor¬ 
dinate  frame  and  then  writing  subroutines  that  will  tabulate  the  coefficients  of  these  differen¬ 
tial  equations  using  the  tabulated  coordinate  transformation.  This  method  has  been  suc¬ 
cessfully  used  in  two  dimensions  so  we  are  concerned  with  the  three  dimensional  extension  of 
the  existing  codes  where  the  major  problem  encountered  is  the  volume  of  algebra  and  coding 
required  to  complete  the  method.  To  overcome  these  difficulties,  a  symbol  manipulation 
program  in  VAXIMA  is  written  that  has  as  input  the  formula  for  the  given  differential 
equation  in  some  natural  coordinates  and  has  as  output  the  required  FORTRAN  subroutines. 
Because  of  the  complexity  of  the  resulting  code,  code  validation  was  performed  by  systematic 
truncation  error  testing.  The  paper  concludes  with  a  discussion  of  the  problems  encountered 
in  using  a  symbol  manipulator  to  write  large  FORTRAN  codes.  ©  >9*5  Academic  Press,  Inc. 

Contents.  1.  Introduction.  2.  Coordinate  transformations.  3.  Difference  equations.  4.  The  sym¬ 
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b.  Numerical  Validation.  6.1.  Hosted  equation  convergence  testing.  6.2.  3D  Hosted  equation 
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The  ELF  (Electric  Field)  codes  have  been 
developed  as  a  design  tool  for  modeling  the 
electric  potential  and  field  strength  in  the 
cavity  of  a  laser  or  switch.  In  an  interactive 
computer  environment,  the  designer  is  able  to 
perturb  the  device  operating  parameters  (applied 
voltage,  external  circuit  elements,  gas  mixtire, 
external  ionization  source  parameters,  etc.) 
and/or  the  electrode  shapes,  and  quickly  obtain 
new  solutions.  The  codes  solve  two-  and  three- 
dimensional  problems  in  both  steady-state  and 
time-dependent  modes.  The  five-year  development 
of  the  codes  has  involved  an  interdisciplinary  team 
approach  of  laser  physicists,  mathematicians,  and 
computational  engineers.  The  computational 
techniques  developed  and  embodied  in  the  codes 
involve  the  following  areas:  semidirect/marching 
methods  for  nonlinear  elliptic  equations,  fully 
implicit  methods  for  strongly  coupled  nonlinear 
time  evolution  equations,  solution-adaptive 
boundary-fitted  grid  generation,  gas  conductivity 
modeling,  parametric  surface  representations, 
super-microcomputer  operations,  artificial 
intelligence  (computer  symbolic  manipulation), 
code  validation  procedures,  software  engineering 
for  interactive  codes,  and  numerical  machining 
considerations,  future  work  may  involve  multi- 
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abstract 


A  single-equation  MG-FAS  solver  is  applied  to  elliptic  gnd  generation  equations 
obtained  from  variational  principles  and  to  problems  of  fluid  dynamics  and  electro¬ 
magnetics  on  these  grids.  The  problem  of  automatic  generation  of  subgnd  coefficients, 
given  the  user  defined  discretized  problem  on  a  defining  grid,  is  addressed  in  a  novel 
and  clear  way  which  has  some  advantages.  The  same  idea  is  applied  to  supergrid 
generation  of  the  elliptic  grid  generation  equations,  which  can  result  in  substantial 
arithmetic  savings  for  solution  adaptive  grid  generation  in  3D. 


1.  INTRODUCTION 

The  motivation  for  the  present  work  on  multigrid  methods  is  the  solution 
of  problems  in  computational  fluid  dynamics  and  electromagnetics,  and  grid 
generation  for  these  problems  [1-4],  in  2D  and  3D.  Several  aspects  of  these 
problems  bear  on  the  formulation  and  the  difficulties  addressed  in  the  present 
work. 

The  well-known  high-grid-Reynolds-number  phenomenon  in  the  CFD 
problems  [5,  6],  and  the  corresponding  phenomenon  arising  from  rapidly 
varying  (or  even  discontinuous)  diffusion  (conductivity)  coefficients  in  plasma 
electromagnetics,  can  give  rise  to  nonphysical  oscillatory  discretized  solutions, 
which  is  aggravated  by  coarse  grid  discretization  within  the  multi  grid  solu¬ 
tion.  The  gnd  generation  equations  can  be  very  complex,  being  based  on  an 
extension  [2,  3]  of  the  Brackbill-Saltzmann  [7]  variational  approach.  We  use 
computer  symbolic  manipulation  to  form  the  variational  equations,  and  to 
write  fortran  subroutines  for  the  stencils  [8.  9].  Even  in  2D.  when  a  strong 
coordinate  angle  control  (nonorthogonal)  is  used  (e.g.  to  adapt  to  solution 
streamline  angles)  the  arithmetic  is  very  expensive,  compared  to  that  of  the 
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SUMMARY 

When  variational  ^nd  generation  techniques  are  used  to  produce  grids  suitable  tor  solving  numerical  partial 
differential  equations  in  irregular  geometries,  a  reference  grid  can  be  used  to  determine  the  properties  ot 
the  desired  grid.  Grid  folding  is  a  major  problem  in  all  methods  of  numerical  grid  generation;  this  paper 
studies  the  use  of  a  reference  grid,  in  variational  smoothing  methods,  to  prevent  such  foldings.  The  effects 
of  various  types  of  reference  grid  are  compared  and  natural  reference  grids  are  shown  to  prevent  foldings  in 
the  e.xamples  studied.  The  analysis  is  performed  for  a  model  problem  with  only  one  free  point,  and  is  then 
apolied  for  finer  grids.  In  addition,  a  short  description  of  a  discrete  variational  method  is  given.  The  reference 
grid  concept  applies  to  this  discrete  formulation. 


INTRODUCTION 

Grid  folding  is  one  of  the  main  concerns  in  numerical  grid  generation.  In  this  paper  we  analyse 
a  variational  grid  generation  method  that  uses  a  reference  grid  as  a  tool  in  preventing  such  folding. 
A  reference  grid  is  a  region  somewhat  like  the  physical  region  where  the  grid  is  to  be  placed  but 
which  can  be  much  simpler  (see  Figure  1).  The  idea  is  to  define,  on  the  reference  grid,  the 
properties  of  the  grid  which  we  want  to  transfer  to  the  given  physical  object.'  (It  is  worth  noting 
that  the  reference  grid  is  not  limited  to  the  determination  of  the  grid  properties  in  the  interior  of 
the  geometric  object,  but  also  can  be  used  to  determine  the  grid  properties  on  the  boundary.') 

In  this  paper  we  are  considering  two  of  the  integrals  from  the  variational  methods  introduced 
bv  oteinberg  and  Roache.'  which  are  similar  to  the  methods  of  Brackbill  and  Saltzman;-  the  two 
functionals  presented  provide  (1)  the  measure  of  spacing  between  the  gnd  lines  (smoothness)  and 
(2)  the  measure  of  the  area  of  the  grid  ceils.  The  minimization  problem  is  usually  solved  by 
calculating  the  Euler-Lagrange  (E-L)  equations  for  the  variational  problem;  the  computer  creates 
a  arid  by  solving  a  centred  finite  difference  approximation  for  the  Euler-Lagrange  equations. 

Our  studv  is  motivated  by  the  fact  that  the  Winslow'  grid  generator  (homogeneous  TTM*) 
produces  folded  grids  in  certain  electrode  configurations.- •  ”  In  these  configurations,  the  regions 
studied  were  roughly  like  the  region  between  two  concentric  ellipses;  the  grid  was  100  <  15.  It 
was  demonstrated  that  the  folding  was  not  due  to  nonlinear  effects  or  coding  error.  The  folding 
is  intrinsic  to  the  Winslow^  (or  homogeneous  T I .VU)  method  when  applied  to  certain  geometries, 
for  finite  grid  sizes.  Although  the  theorem  of  .Mastin  proves  that  the  grids  will  not  fold  in  two 
dimensions'^  or  three  dimensions'*  in  the  continuum  limit,  the  results  of  Roache  and  Steinberg- 
conclusively  and  succinctly  demonstrated  that  the  theorem  does  not  apply  for  finite  mesh 
increments.  The  difficulty  manifested  itself  even  in  the  simplest  possible  grids  and  the  results  for 
the  simple  grids  proved  relevant  to  more  complex  and  finer  resolution  cases  (see  Reierences  5.  9 
and  10).  When  only  3x3  gnds  are  considered,  as  the  case  for  the  present  analyses,  the  region 
degenerates  to  the  area  between  two  triangles,  as  in  Figure  3. 
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Abstract 

In  numerical  grid  generation,  having  a  grid  fold  is  a  major  problem.  This  paper  studies 
techniques  for  preventing  such  grid  foldings.  The  analysis  is  done  for  a  model  problem  with 
the  simplest  possible  numerical  grid,  one  where  there  is  only  one  point  to  be  determined. 

1  Introduction 

Our  main  interest  is  to  study  the  grid  generation  techniques  introduced  by  Steinberg  and 
Roache  in  [l,  2|.  These  techniques  are  variational  and  are  closely  related  to  the  techniques 
introduced  by  Brackbill  and  Saltzman  in  [3].  We  will  also  compare  the  variational  methods 
to  a  simple  version  of  a  method  introduced  by  Thompson,  et.al  [4].  There  are  three  basic 
properties  to  be  controlled  in  a  grid:  smoothness,  cell  area  or  volume,  and  orthogonality. 
These  properties  are  represented  by  integrals  which  are  to  be  minimized.  The  derivation 
of  these  integrals  has  been  studied  in  several  places  [l,  3’. 

Our  study  is  motivated  by  the  discovery  that  the  homogeneous  TTM  or  Winslow  grid 
generator  produced  folded  grids  in  certain  electrode  configurations  [2’.  In  these  configu¬ 
rations  the  regions  studied  were  roughly  like  the  region  between  two  concentric  ellipses 
[see  Figure  laj.  In  the  electrode  problem  the  folded  grid  was  100  by  15.  It  was  discovered 
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SUMMARY 

A  hybrid  technique  for  adaptive  grid  generation  using  a  variety  of  Poisson  grid  generators  is  described. 
The  technique  ensures  that,  when  adaptivity  functions  are  weak,  the  adapted  grid  reverts  to  an  arbitrary 
specified  base  grid  rather  than  to  the  grid  produced  by  the  homogeneous  elliptic  generator.  The  hybrid 
technique  also  serves  to  e.xpose  the  weakness  of  the  claim  that  elliptic  grid  generators  produce  smooth 
grids. 


INTRODUCTION 

Adaptive  grid  generation  has  been  used  frequently  in  recent  years  by  investigators  to  obtain 
more  accurate  and  efficient  solutions  to  a  variety  of  problems.  This  paper  identifies  one 
common  shortcoming  of  adaptive  elliptic  methods,  and  presents  a  hybrid  technique  to 
overcome  it.  This  simple  technique  also  explains  a  common  misconception  regarding  grid 
smoothness. 


ADAPTIVE  POISSON  GRID  GENERATION 

Using  the  idea  of  equidistribution  of  some  weight  functions  over  the  computational  mesh, 
Anderson’  showed  that  the  Poisson  grid  generation  scheme  of  Thompson,  Thames  and 
Mastin'  (denoted  as  TTM  herein)  can  be  interpreted  as  a  non-linear  equidistribution  scheme 
and  extended  TTM  to  include  adaptivity.  A  variety  of  grid  controls  may  be  accommodated, 
with  the  grid  control  functions  denned  in  terms  of  auxiliary  functions  introduced  by  Thomas 
and  Middlecoff.  ^  The  method  allows  a  relatively  simple  conversion  of  the  ubiquitous  TTM 
codes  to  solution  adaptivity  through  the  non-homogeneous  terms. 

However,  this  adaptive  method  (and  others  based  on  elliptic  grid  generators)  has  a  serious 
shortcoming.  Whenever  the  weighting  function  is  weak,  the  generated  grid  reduces  to  that 
produced  by  the  original  WinsioW*  or  homogeneous  TTM'  method.  Thus,  the  grid  will 
e.xperience  the  same  type  of  difficulties  as  that  method.  Thompson  ero/.‘  noted  the  tendency 
of  the  homogeneous  method  to  pull  grid  lines  away  from  any  concave  boundary.  Even  more 
seriously,  Roache  and  Steinberg  demonstrated  that,  for  a  pamcuiar  common  ciass  of 
geomet.'-y,  the  homogeneous  method  will  produce  folded  grids.  (See  also  the  e,xample 
calculations  in  the  original  Winslow*  paper.) 


*  Also,  Professor,  Department  of  .Mathematics  and  Statistics.  University  of  New  Mexico. 
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Variational  Grid  Generation' 
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Recently,  variational  methods  have  been  used  to  numencally  generate  grids  on  geomet- 
ometric  objects  such  as  plane  regions,  volumes,  and  surfaces.  This  amcte  presents  a 
new  method  of  determining  variational  problems  that  can  be  used  to  control  sucn  proper¬ 
ties  of  the  2nd  as  the  spacing  of  the  points,  area  or  volume- of  the  cells,  and  the  angles 
between  the  gnd  lines.  The  methods  are  applied  to  curves,  surfaces,  and  volumes  in 
three-dimensional  space:  then  segments,  plane  curves,  and  plane  regions  appear  as  spe¬ 
cial  cases  of  the  general  discussion.  The  methods  used  here  are  simpler  and  clearer  and 
provide  more  direct  control  over  the  grid  than  methods  that  appear  elsewhere.  The 
methods  are  applicable  to  any  simply  connected  region  or  any  region  that  can  be  made 
simply  connected  by  inserting  artificial  boundaries.  The  methods  also  generalize  easily 
to  solution-adapuve  methods. 

An  important  ingredient  in  our  method  is  the  notion  of  a  reference  grid.  A  reference 
grid  is  defined  on 'a  region  that  is  simpler,  but  analogous  to.  the  geometric  object  on 
which  a  grid  is  desired.  Variational  methods  are  then  used  to  transfer  the  reference  gnd 
to  the  geometric  object.  This  gives  simple  and  precise  control  of  the  local  properties  of 
the  grid. 


I.  INTRODUCTION 

In  this  article  we  study  the  problem  of  using  a  computer  program  to  automati¬ 
cally  generate  grids  on  curves,  on  surfaces,  or  in  volumes  in  three-dimensional 
Euclidean  space.  Curves  and  surfaces  in  two-dimensional  Euclidean  space  and 
segments  will  be  included  as  special  cases  in  the  general  discussion.  The  prob¬ 
lem  of  grid  generation  has  a  long  history  and  many  important  applications, 
which  are  described  in  [4],  We  were  motivated  to  smdy  the  grid  generation 
problem  because  we  were  interested  in  using  such  ends  in  finite-difference 
codes  that  are  used  to  solve  partial  differential  equations.  However,  that  appli¬ 
cations  area  does  not  play  a  central  role  here. 

Intuidveiy,  a  discrete  grid  on  a  curve  is  a  sequence  of  points  that  divides  the 
curve  into  pieces  that  are  nearly  straight  line  segments;  a  discrete  gnd  on  a  sur¬ 
face  is  a  set  of  points  that  divides  ±e  sunace  into  approximate  parallelograms 
(the  grid  need  not  be  orthogonal);  a  discrete  gnd  in  a  volume  is  a  set  of  pomts 
that  divides  the  volume  into  ceils  that  are  nearly  parallelepipeds .  In  all  cases. 


^Xhis  work  was  partially  supported  by  the  U.S.  Air  Force  Office  of  Scientific  Re¬ 
search.  by  the  U.S.  Anny  Research  Office,  and  by  the  National  Science  Foundation 
Grant  MCS-9102683. 
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Abstract 

Two  approaches  to  numerical  variational  grid  generation  ’n-e  a 
variational  problems  to  control  the  properties  of  the  grid.  ?r 
the  linear  combination  were  determined  by  experimentation, 
parameters  can  be  obtained  from  simple  model  problems. 

*Tliis  work  was  partially  supported  by  ARO  and  AFOSR. 
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.Abstract.  The  symbol  manipulator  MACSYMA  is  used  as  a  high-level  interface  between 
the  user  and  FORTRAN.  MACSYMA  is  used  to  write  numerical  code  for  solving  sys¬ 
tems  of  ordinary  differential  equations. 

Kes-ivords.  Symbol  manipulation:  automatic  code  generation. 


INTRODUCTION 

Languages  such  as  FORTRA.N  provide  valuable  tools 
to  the  computer  user,  and  such  •‘mid-level"  languages 
can  be  used  to  do  sophisticated  .numerical  computing. 
However,  there  are  some  problems  which  arise.  One  is 
that  in  FORTR.AN  all  code  is  specific  to  the  calculation 
at  hand:  array  sizes  must  be  set  up  in  the  code,  and 
input  to  pre-written  codes  is  severely  limited.  With 
FORTR.A.V  code,  it  is  basically  ■‘numbers  in.  numbers 
out  ’.  Equations  cannot  be  changed  without  editing  the 
code.  This  means  that  to  write  codes  for  solving  many 
different  systems  of  the  same  type,  the  user  must  maxe 
a  -tempiate  file'’,  copy  it,  and  insert  ail  array  sizes  and 
equations  individually,  w'nich  is  highly  inefficient.  .An¬ 
other  problem  is  that  input  and  output  for  rORTR.AN 
codes  are  .not  easiiy  interpreted  by  those  not  familiar 
•'vit.h  the  caicuiations  in  question.  The  user  must  often 
process  output  to  make  it  understandable  to  others. 

Both  of  these  problems  can  be  eliminated  by  the  'use 
of  a  high-level  interface  between  users  and  FORTRAN 
code.  The  plan  here  is  for  the  user  to  type  in  readily- 
understood  expressions  -  statements  and  equations  - 
which  are  used  by  a  symbol  mamipuiator  ;e.g.  M.AC- 
SN'NlAl  to  generate  FORTR.AN  code:  this  in  turn  can 
be  used  to  drive  some  standard  FORTR.A.N  code.  Much 
of  the  intermediate  information  .MACSYMA  produces 
on  its  way  to  the  FORTRAN  code  can  be  ontpnt  to  the 
user,  who  can  then  see  resnlts  in  the  form  of  eqnations. 
rather  than  lists  of  numbers.  Output  from  FORTRAN 


code  can  also  be  routed  directly  to  graphics  codes,  to 
produce  plots  which  give  a  much  clearer  picture  of  the 
results  than  lists  of  numbers  could.  The  ultimate  goal 
is  to  begin  with  ■‘high-level”  mathematics,  physics,  and 
engineering,  and  produce  results  that  are  easily  under¬ 
stood  by  humans. 

The  purpose  of  this  paper  is  to  show,  using  a  simple 
example,  .tow  M.ACSYNIA  can  'oe  used  as  an  interface 
between  the  user  and  ■‘canned"  FORTR.AN  codes  lin 
t.Pis  case  t.ne  code  RKr45,  a  numerical  ODE  integra¬ 
tor'!  . 

PRELIMIN  ARIES 

Before  ioo'xing  at  -..he  MACSYNLA  code,  it  is  heipfui  to 
know  a  .ittle  about  a  few  MACSYNt.A  commands,  and 
.now  NLACS'AR.I.A  is  programmed.  This  section  gives  a 
i^rief  expianation  of  some  of  the  commands  used  in  the 
codes  whic.h  follow. 

One  of  the  data  types  available  in  MACS'VNIA  is  the 
list.  One  way  to  make  a  list  is  to  use  the  command 
tnacoTu.  This  J  done  by  producing  an  empty  iist.  and 
then  adding  elements  to  the  end  of  the  iist: 

Cc_l)  lisc  :  []  ; 

II  («i.) 


1107 


J.  Symbolic  Computation  (1986)  2,  213-216 


Using  MACSYMA  to  Write  FORTRAN  Subroutinest 

STANLY  STEINBERG 

Department  of  Mathematics  and  Statistics. 

University  of  .\'e\v  Mexico.  .Albuquerque  .VA/87131. 

PATRICK  J.  ROACHE 

EcoJ\  namics  Researci:  Assoctates.  Inc.. 

P.O.  Box  SI 72,  .Albuqucraue  SM  S"198.  U.S..A. 

(Received  2  January  1986) 


1.  Introduction 

This  letter  describes  some  ongoing  work  that  uses  the  symbol  manipulator  macsyma  to 
wnte  FORTRAN  subroutines  which  are  incorporated  into  a  large  finite  difference  code  that 
is  used  to  solve  boundary  value  problems  for  elliptic  partial  differential  equations.  The 
boundary  value  problems  model  the  steady  state  of  certain  physical  devices  such  as  the 
electric  field  in  a  laser  cavity  or  the  flow  around  a  blade  in  a  turbine.  The  devices  that 
interest  us  are  those  that  occupy  a  fairly  complicated  region  in  three-dimensional  space 
or.  at  best,  have  sufficient  symmetry  so  that  they  can  be  modelled  using  a  complicated 
two-dimensional  region.  The  partial  differential  equations  that  are  used  to  model  the 
physics  of  the  device  are  called  the  hosted  equations.  The  hosted  equations  tend  to  be 
rather  simple  non-linear  elliptic  equations.  However,  our  methods  will  apply  to  the  most 
complicated  hosted  equations.  The  problem  is  complete  when  appropriate  boundary 
conditions  are  specified.  These  problems  are  difficult  to  solve  because  of  the  complicated 
geometry  of  the  region. 


2.  Overview  of  the  Problem 

Several  important  approaches  to  solving  such  problems  are  based  on  finding  a  change 
of  coordinates  that  maps  the  given  physical  region  into  a  rectangular  region  in  logical 
space.  Such  a  coordinate  change  then  generates  a  grid  in  the  physical  region  that 
corresponds  to  a  rectangular  grid  in  logical  space.  This  approach  now  has  a  long  history 
that  IS  referenced  in  the  papers  listed  below.  Once  such  a  transformation  has  been  found, 
the  hosted  equation  and  the  boundary  conditions  are  transformed  to  the  new  coordinate 
system  and  then  standard  finite  difference  methods  are  applied  to  solve  the  problem. 
These  numencal  methods  are  described  in  our  papers.  What  we  will  descnbe  here  is  how 
a  symbol  manipulator  is  used  to  generate  the  fortran  subroutines  that  are  used  by  the 
finite  difference  codes. 

The  regions  that  interest  us  are  so  complicated  that  it  is  often  impractical  or  impossible 
to  find  an  analytic  change  of  coordinates  that  transforms  the  region  to  a  rectangular 
region.  Historically,  practitioners  in  this  area  have  used  finite  difference  approximations 

+  Work  supponed  partially  by  the  U.S.  Army  Research  Office  and  the  U.S.  .Air  Force  Office  of  Saentific 
Research. 
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In  addition  to  the  publications  listed  previously, 
which  were  directly  supported  by  the  subject  contract,  the 
continuation  of  this  grid  generation  work  initially  funded 
by  AFOSR  has  resulted  in  other  significant  publications. 
Besides  the  two  books  mentioned  earlier,  a  sampling  of 
related  later  work  follows,  especially  those  resulting  from 
the  AFWL  contracts  entitled  "Adaptive  Grids  for  RAVEN"  (see 
Section  G  below) . 
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Attempts  to  use  vanational  grid-generation  methods  to  generate  gnds  on  certain  surfaces 
of  modest  shape  failed.  There  were  suilicient  points  in  the  ends  to  well-resolve  the  surface,  so 
the  failures  were  not  easily  explained.  Similar  dilTiculties  were  found  for  variational  end 
generation  on  cuiwes;  those  problems  are  caused  by  multiple  solutions  of  the  underlying  non¬ 
linear  algebraic  equations.  C  I9«0  .Xcaaemic  Press.  Inc. 


1.  Backgrount) 

This  paper  describes  several  grid-generation  anomalies  discovered  by  the  authors 
while  using  the  vanational  techniques,  descnbed  in  Steinberg  and  Roache  [9],  to 
generate  grids  on  surfaces.  These  surfaces  are  of  modest  shape  (they  are  intended 
to  model  a  stem  wave  behind  a  ship  hull)  and  the  grids  being  generated  resolve  the 
surfaces  well.  The  anomalies  manifest  themselves  as  convergence  problems.  Grids 
are  difficult  or  impossible  to  generate;  for  some  coarse  resolutions  ito  be  e.xpected); 
for  many  fine  resolutions  mot  to  be  e.xpected).  For  many  intermediate  resolutions, 
it  is  easy  to  generate  excellent  grids.  Numencai  expenmentation  fails  to  pinpoint 
the  source  of  the  difficulties.  The  surface  grid-generation  equations  are  a  com¬ 
plicated  system  of  two  coupled  quasi-linear  partial  differential  equations  that  have 
a  form  similar  to  elliptic  equations. 

Some  numerical  experimentation  shows  that  the  analogous  gnd-generation  code 
for  curves  has  the  same  diflicuities  as  the  surface  gnd  generator.  This  paper  presents 
an  analysis  of  the  curve  problem  because  it  is  simpier  than  the  surface  problem  and 


•  This  work  was  partially  supported  by  the  Office  of  Naval  Research  and  the  .Air  Force  Weapons 
Laboratory.  This  is  a  modified  version  of  AIAA  Paper  No,  ,S8-3T40.  presented  at  the  1st  National  Ruid 
Dynamics  Conference.  Cinannati.  OH.  July.  1988. 
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Need  for  Control  of  Numerical  Accuracy 


Patrick  J.  Roache* 

Ecodynamics  Research  Associates,  Albuquerque,  New  Mexico  87198 


The  need  for  the  control  of  numerical  accuracv  in  computational  fluid  dynamics  (CFD)  code  solutions  is 
reviewed  in  the  light  of  current  journal  practice  and  experience  with  implementation  of  an  editorial  policy  on  the 
same  subject  published  by  the  Journal  of  Fluids  Engineering.  Various  actual  objections  to  that  policy  are  listed 
and  responses  are  given.  The  general  successes  and  particular  difficulties  experienced  in  the  implementation  of 
the  policy  are  noted.  The  broader  question  of  code  verification,  validation,  and  certification  is  considered.  It  is 
suggested  that  professional  societies  such  as  the  .Vl.4.4  and  .American  Society  of  Mechanical  Engineers  may 
ultimately  become  involved  in  the  task  of  certification  of  commercially  available  CFD  codes. 


I.  Introduction 

HIS  paper,  on  the  general  philosophy  and  need  for 
numerical  accuracy  control  in  computational  fluid  dy¬ 
namics  (CFD)  codes,  is  based  on  experience  with  the  imple¬ 
mentation  of  an  editorial  policy  statement  by  the  American 
Society  of  Mechanical  Engineers  (ASME),  published  in  the 
Journal  of  Fluids  Engineering  iJFEiJ  The  policy  statement 
was  conceived  following  the  creation  of  the  position  of  .Associ¬ 
ate  Editor  for  Numerical  .Methods  in  the  JFE.  which  formally 
recognized  the  special  needs  of  this  discipline.  The  .ASME 
policy  and  supporting  statements  are  reproduced  in  the  .Ap¬ 
pendix.  The  rationale  and  needs  for  the  policy  statement  are 
explained  in  the  announcement.  The  JFE  had  previously  pub¬ 
lished  and  had  many  years  of  experience  with  a  similar  re¬ 
quirement  for  uncertainty  analysis  in  experimental  papers. 
Following  our  early  e.xperience  with  this  new  policy,  a  similar 
policy  was  also  adopted  by  the  ASME's  Journal  of  Heat 
Transfer.  The  general  subject  of  control  of  numerical  accu¬ 
racy  has  become  something  of  a  ‘‘hot  topic."  with  special 
reference  to  the  National  Aero-Space  Plane,  a  session  at  the 
AIA.A  1989  Thermophysics  Conference.  .ASME  sessions  at  the 
1988  and  1989  Winter  .Annual  Meetings,  and  by  the  Texas 
Institute  for  Computational  .Mechanics  workshop  on  the 
slightly  broader  topic  of  reliability  in  computational  mechan¬ 
ics  m  October  1989. 

II.  Resistance  and  Objections 

Although  Roache  et  al.'  though:  that  the  policy  statement 
and  the  editorial  requirement  were  quite  mild,  it  was  not 
universally  welcomed.  Objections  were  offered  by  some  of  the 
other  editorial  board  memoers  and  bv  other  members  of  the 
professional  community  from  whom  1  solicited  comments  in 
the  months  following  the  publication  of  the  statement.  Some 
of  these  objections,  all  of  w  hich  are  actual  li.e..  not  straw-man 
objections),  are  listed  below,  together  with  my  responses. 

Objection  t 

It  IS  too  expensive  of  computer  time  to  do  mesh  doubling 
calculations  in  order  to  ascertain  grid  convergence. 

There  are.  of  course,  other  ways  to  ascertain  grid  conver¬ 
gence  than  the  straightforward  method  of  grid  doubling.  In 


Presented  as  Paper  89-1669  at  the  24th  AlAA  Thermophysics  Con¬ 
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fact,  this  is  probably  the  most  reliable  method  available,  but 
there  are  other  approaches,  as  briefly  touched  upon  in  the 
original  policy  statement.  If  the  cost  of  computer  resources  is 
not  a  problem  to  the  researcher,  this  is  certainly  the  easiest 
approach  to  take.  However,  if  computer  resources  are  a  prob¬ 
lem.  there  are  other  methods  that  are  not  intensive  users  of 
computer  time.  It  seems  that  the  greater  objection  for  doing 
grid  convergence  studies  is  the  fact  that  it  requires  a  bit  of 
conscientious  work  on  the  pan  of  researchers. 

Also,  if  it  is  argued  that  it  is  simply  too  expensive  to  do  any 
kind  of  control  of  numerical  accuracy,  then  1  would  argue  that 
the  author  simply  cannot  be  in  this  CFD  business.  .After  all.  if 
you  do  not  have  a  wind  tunnel  you  cannot  do  experimental 
testing.  My  impression  of  the  situation  is  actually  worse  than 
this.  Journal  articles  in  the  late  1960s  and  early  1970s  com¬ 
monly  predicted  high  resolution  accuracy  runs  when  the  next 
generation  of  computers  became  available.  But  for  the  most 
pan  they  have  not  been  used  that  way.  Generally  the  tremen¬ 
dous  advance  of  computing  power  has  been  used  to  produce 
more  mediocre  papers  rather  than  fewer  reliable  ones. 

Objection  2 

Some  e.vception  to  the  policy  should  be  made  for  e.vpensive 
calculations,  particularly  three-dimensional  turbulent  studies. 

1  do  not  think  that  the  overall  cost  of  the  computations 
should  be  a  consideration.  It  seems  clear  tnat  the  incremental 
cost  of  performing  a  grid  convergence  test  should  be  normal¬ 
ized  by  the  cost  of  the  base  case.  In  this  sense  it  is  cheaper  to 
validate  the  grid  convergence  on  a  three-dimensional  problem 
than  on  a  two-dimensional  probietn,  presuming  that  the  incre¬ 
mental  work  involves  an  extra  coarse  grid  computation.  This 
again  relates  to  the  first  point,  which  states  tnat  ;;  is  not 
necessary  to  do  a  grid  doubling  in  order  to  ascertain  some 
index  of  numerical  accuracy.  A  grid  halving  is  also  aporopri- 
atc.  Of  course,  the  advantage  lies  in  doing  a  end  doubiing  test 
because  the  error  bounds  will  be  sharper. 

Objection  3 

Turbulence  modeling,  rather  than  the  numerical  solution  of 
the  partial  differential  equations,  is  the  real  determinator  of 
accuracy. 

My  response  is,  yes  and  no.  Accuracy  is  a  question  to  be 
addressed  even  for  laminar  flow  calculations,  wherein  the 
constitutive  equations  are  not  in  doubt.  The  discretization 
error  does  not  disappear  just  because  one  uses  a  turbulence 
model!  Our  point  in  the  JFE  policy  statement,  and  the  First 
criticism  which  our  evaluation  committee  made  at  the  1980/81 
Stanford  meeting  on  comple.x  turbulence  flows,-  is  that  one 
cannot  evaluate  different  turbulence  models  unless  one  first 
satisfies  grid  convergence.  There  are  yet  more  considerations 
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Abstract 


The  influence  of  sweep  on  deep  dynamic  stall  of  a  rapidly 
pitching  swept  wing  at  low  Mach  number  with  laminar  flow  has  been 
investigated  through  the  use  of  numerical  flow  simulations.  The 
problem  involves  the  modeling  of  a  wind  tunnel  test  section  in 
which  the  wing  spans  the  tunnel.  The  flow  Reynolds  number  is 
10,000,  free  stream  Mach  nximber  is  0.2,  the  reduced  frequency  is 
equal  to  0.3,  and  the  sweep  angle  is  30  degrees.  The  solution  of 
the  full  unsteady  three-dimensional  compressible  Navier  Stokes 
equations  was  obtained  on  the  CRAY-2  supercomputer  through  use  of 
an  implicit  finite  difference  Approximate  Factorization  algorithm 
coupled  with  a  non-orthogonal  moving  grid.  The  sweep  effects 
have  been  determined  by  comparing  the  unswept  (2-D)  and  swept 
(3-D)  solutions.  Sweep  tends  to  delay  the  onset  of  dynamic  stall 
and  reduce  the  magnitude  of  unsteady  aerodynamic  loads.  However, 
the  intensity  of  these  effects  vary  significantly  along  the  span 
of  the  wing.  Sweep  was  also  found  to  inhibit  the  growth  of  the 
secondary  stall  vortex  and  its  contribution  to  aerodynamic 
forces.  The  side  tunnel  walls  tend  to  alter  the  behavior  of  the 
stall  vortex  or  prevent  its  formation. 


I .  Introduction 


History 

Around  1963  design  engineers  in  the  helicopter  industry 
obser/ed  that  the  lift  generated  by  a  helicopter  blade  in 
experiments  was  much  higher  than  the  lift  predicted  using 
conventional  aerodynamics.  Harris  and  Pruyn  [1]  realized  the 
importance  of  unsteady  aerodynamics  in  resolving  this 
discrepancy.  Concurrently,  Ham  and  Carelick  [2]  observed  that 
the  extra  lift  was  associated  with  a  vortex  which  formed  by  rapid 
pitching  of  an  airfoil  during  the  unsteady  motion.  Ham  [3] 
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Most  3D  finite  difference  or  rintte  element  codes  make  use  ot  hexahedral  ceils  having  six  faces,  eight 
corners  and  twelve  edges.  Interpolation  between  grids  makes  use  of  the  isoparametric  map  from  a 
uniform  logical  space  to  a  given  hexahedral  cell.  In  two  dimensions,  it  is  well-known  that  the 
isoparametric  map  is  invertible  if  and  only  if  the  Jacobians  at  the  four  cell  corners  are  positive.  The 
corresponding  conjecture  in  three  dimensions  is  shown  to  be  false  and  alternative  tests  for  positivity  of 
the  Jacobian  are  investigated. 


I.  Properties  of  the  Jacobian  in  K' 

Following  Strang  and  Fix  [1],  the  isoparametric  mapping  from  the  unit  square  L\_  = 
{{£.  r})\0^  Tj  ^  1)  to  a  quadrilateral  Q  in  R‘  (Fig.  1)  is  given  by  the  pair  of  bilinear  forms 

v)  =  -^I  +  {■'^2  -  +  (^3  “  (-^-I  ~  •  (1) 

v(  v)  =  >’i  +  (v;  -  y,)f  +  (y-  -  y,)T7  +  (y,  -  y-  -  y,  ^  .  (2) 

The  isoparametric  map  is  a  special  case  of  the  well-known  ’shearing  transformation"  used  in 
algebraic  grid  generation  in  which  all  four  sides  of  the  given  region  are  straight  lines.  The  grid 
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Abstract 


Abstract:  The  Richardson  extrapolation  method,  which 

produces  a  4th  order  accurate  solution  on  a  subgrid  by  combining 
2nd  order  solutions  on  the  fine  grid  and  the  subgrid,  is 
"completed",  in  the  sense  that  a  higher  order  accurate  solution  is 
produced  on  all  the  fine  grid  points. 


INTRODUCTION 

In  his  classic  paper  in  1910,  Richardson  [1]  presented  a 
method  for  obtaining  4th  order  accurate  solutions.  The  method, 
known  variously  as  Richardson  extrapolation,  extrapolation  to  the 
limit,  deferred  approach  to  the  limit,  or  iterated  extrapolation 
[2]  takes  separate  2nd  order  solutions  on  a  fine  grid  and  on  the 
subgrid  formed  of  alternate  points,  and  combines  them  to  obtain  a 
4th  order  solution  on  the  subgrid.  It  is  also  the  basis  of 
Romberg  integration  [ 3 ] . 

The  usual  ass;imptions  of  smoothness  apply,  as  well  as  the 
assumption  (or  perhaps  presumption)  common  to  finite  difference 
methods  that  the  local  error  is  indicative  of  global  error.  The 
method  must  be  used  with  considerable  caution,  since  it  involves 
additional  assumptions  of  monotone  truncation  error  convergence  in 
the  mesh  spacing  h  (which  may  not  be  valid  for  coarse  grids)  and 
since  it  magnifies  machine  round-off  errors  and  incomplete 
iteration  errors  [2,4].  In  spite  of  these  caveats,  the  method  is 
extremely  convenient  to  use  compared  to  forming  and  solving  direct 
4th  order  discretizations,  which  involve  more  complicated 
stencils,  wider  band-width  matrices,  special  considerations  for 
near-boundary  points  and  non-Dirichlet  boundary  conditions, 
additional  stability  analyses,  etc.,  especially  in  nonorthogonal 
coordinates  which  generate  cross-derivative  terms  and  generally 
complicated  equations.  Such  an  application  was  given  in  [5]  by 
the  first  author.  The  method  is  in  fact  oblivious  to  the 
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E.  PROFESSIONAL  PERSONNEL  ASSOCIATED  WITH  THE  RESEARCH 
EFFORT 


Dr.  Patrick  J.  Roache,  Principal  Investigator 
Dr.  Stanly  Steinberg 
Dr.  Jose'  Castillo 
Mr.  Paul  Sery 

F.  INTERACTIONS  (COUPLING  ACTIVITIES) 

(1)  PRESENTATIONS  AT  MEETINGS,  CONFERENCES,  SEMINARS 


PRESENTATIONS 


"Computers  vs.  Algorithms",  Symposium  on  the  Impact  of  Large-Scale 
Computing  on  Air  Force  Research  and  Development,  4-6  April  1984, 
Air  Force  Weapons  Laboratory,  Albuquerque,  NM. 

"Interactive  Electric  Field  Calculations  for  Lasers",  AIAA  17th 
Fluid  Dynamics,  Plasma  Physics,  and  Lasers  Conference,  25-27  June 
1984,  Snowmass,  CO. 

"Computers  vs.  Algorithms",  Engineering  Faculty  Seminar,  5  November 
1984,  University  of  Canterbury,  Christchurch,  New  Zealand. 

"Symbolic  Manipulation  and  Computational  Fluid  Dynamics", 
Engineering  Faculty  Seminar,  7  November  1984,  University  of 
Canterbury,  Christchurch,  New  Zealand. 

"Application  of  a  Single-Equation  MG-FAS  Solver  to  Elliptic  Grid 
Generation  Equations  (Sub-grid  and  Super-grid  Coefficient 
Generation)",  Second  Copper  Mountain  Conference  on  Multigrid 
Methods,  1-3  April  1985,  Copper  Mountain,  CO. 

"Symbolic  Manipulation  and  Computational  Fluid  Dynamics", 
Engineering  Faculty  Seminar,  2  May  1985,  Univ.  California  at  Davis. 

"Symbolic  Manipulation  and  Computational  Fluid  Dynamics",  Applied 
Mathematics  Seminar,  7  May  1985,  Sandia  Labs. ,  Livermore,  CA. 

"A  New  Approach  to  Grid  Generation  Using  a  Variational 
Formulation",  AIAA  7th  Computational  Fluid  Dynamics  Conference,  15- 
17  July  1985,  Cincinnati,  OH. 

"The  ELF  Codes;  Electrode  Design  for  Lasers  and  Switches",  Invited 
Paper,  CTAC-85  Conference,  25-28  August  1985,  Melbourne.  Australia. 

"PDE  Discretization  on  Subgrids  and  Supergrids",  SIGNUM  Meeting,  8 
December  1985,  Albuquerque,  NM. 

"A  Tool  Kit  of  Symbolic  Manipulation  Programs  for  Variational  Grid 
Generation",  AIAA  Aerospace  Sciences  Meeting,  6-9  January  1986, 
Reno,  NV. 

"Symbolic  Manipulation  and  Computational  Fluid  Dynamics",  Invited 
Presentation,  NSF  Workshop  on  Computational  Engineering,  24-25  June 
1986,  NSF  San  Diego  Supercomputer  Center,  San  Diego,  CA. 

"The  ELF  Codes",  Spectra  Technologies,  Seattle,  WA,  17  July  1986. 

"Grid  Generation  in  Finite  Analysis",  Applied  Mechanics  Research 
Seminar,  Mechanical  Engineering  Department,  University  of  New 
Mexico,  February  17,  1987,  Albuquerque,  NM. 

"Symbolic  Generation  of  Strong  Conservation  Forms  for  General 
Coordinate  Finite  Difference  Codes",  Invited  Presentation, 
Symposium  on  Symbolic  Computation  in  Heat  Transfer  and  Fluid 
Mechanics,  ASME  Winter  Annual  Meeting,  Chicago,  Ill.  Dec.  1988. 
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(1)  CONSULTATIVE  AND  ADVISORY  FUNCTIONS 

P.  J.  Roache  and  S.  Steinberg  consulted  with  Dr.  R. 
L.  Rapagnani  and  other  personnel  at  AFWL  on  various 
problems  involving  grid  generation  during  1984-1986. 

G.  NEW  DISCOVERIES,  INVENTIONS,  PATENTS,  SPECIFIC 
APPLICATIONS 

No  inventions  or  patents  have  resulted  from  this  work. 

The  2-D  and  3-D  adaptive  grid  generation  algorithms 
developed  in  this  work  and  the  preceding  AFOSR  contract  have 
been  applied  in  the  ELF  codes  (ELectric  Field)  developed  at 
Ecodynamics  under  subcontract  to  Tetra  Corporation  for  the 
Air  Force  Weapons  Laboratory  (now  Phillips  Laboratory, 
Kirtland  AFB,  NM) .  These  codes  are  under  fairly  extensive 
use  at  Phillips,  at  Tetra  Corporation,  and  at  university  and 
for-profit  DoD  contractors  for  the  design  and  optimization 
of  laser  electrodes  and  high  power  switches.  Applications 
are  in  the  Pulsed  Power  area,  including  SDI  research  and 
development.  Over  25  copies  of  ELF  have  been  sold 
commercially,  many  to  government  and  Air  Force  contractors, 
as  well  as  the  Canadian  Atomic  Energy  Commission.  The  ELF 
codes  were  also  the  utilized  in  an  SBIR  Phase  I  contract 
from  AFOSR  to  Ecodynamics  on  "Design  Optimization  of  Systems 
Governed  by  Partial  Differential  Equations". 

The  2-D  and  3-D  adaptive  and  variational  grid 
generation  algorithms  developed  in  this  work  were  also  the 
basis  for  SBIR  Phase  I  and  II  contracts  to  Ecodynamics  from 
the  Air  Force  Weapons  Laboratory  (now  Phillips  Laboratory, 
Kirtland  AFB,  NM)  entitled  "Adaptive  Grids  for  RAVEN".  The 
IMPLICIT  RAVEN  code  developed  under  this  contract,  as  well 
as  the  earlier  RAVEN  code,  are  accessible  on  the  Cray-2 
system  library  at  Phillips  and  are  available  for  use  by 
government  contractors  for  the  analysis  of  problems  in 
subsonic  and  supersonic  aerodynamics  and  well  as  chemical 
lasers.  The  IMPLICIT  RAVEN  code  has  also  been  marketed 
commercially. 

These  grid  generation  algorithms  were  also  the  basis 
for  SBIR  Phase  I  and  II  contracts  to  Ecodynamics  from  the 
U.S.  Army  (AMSAV)  on  "Dynamic  Stall  Control",  specifically 
focused  on  retreating  helicopter  blade  dynamic  stall. 

(Phase  II  has  just  recently  started.) 

They  were  also  the  basis  for  a  contract  from  the  U.S. 
Navy  (ONR)  on  "Computational  Ship  Waves"  with  specific  focus 
on  free  surface  calculations  behind  a  transom  stern.  Dr. 

R.  Yeung  of  University  of  California  at  Berkeley  and  his 
graduate  students  have  highly  tuned  our  variational  surface 
grid  generation  method  (utilizing  a  reference  grid  to 
control  target  grid  properties)  for  free  surface  flows  with 
great  success,  and  have  published  their  results  recently  in 
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International  Journal  for  Numerical  methods  in  Fluids. 

These  grid  generation  algorithms  were  also  the  basis 
for  Phase  I  SBIR  contracts  from  the  U.S.  Navy  (NUSC)  on 
•'Unsteady  Flow  in  a  Torpedo  Shuttleway",  from  the  U.S.  Navy 
(NAVSEA)  on  "Computational  Modeling  of  3-D  Unsteady  Flow", 
specifically  focused  on  simulation  of  flow  in  gas  turbine 
diffusers,  and  from  the  National  Science  Foundation 
(Mathematics  Division)  on  "Robust  and  Fast  Numerical  Grid 
Generation".  The  Phase  II  award  decisions  for  the  last  two 
have  not  yet  been  made. 

Ecodynamics  is  presently  planning  to  market  these  grid 
generation  modules  interfaced  to  the  two  most  used  packages 
in  the  aerodynamics  community,  GRIDGEN  and  EAGLE. 


